%**************************************************************************
%
%           "Endogenous Liquidity and Capital Reallocation"
%      Cui, W., Wright, R., & Zhu, Y. (2024), Journal of Political Economy
%                       Last Modified: Feb 2024.  
%
%**************************************************************************

%%%% Computing steady state and comparative statics with search frictions

clc; clear; close all; format short

%% 1. --------------------      Initilization    --------------------------

%%%% This step is important to initiate if parameters were not saved
%Main_iid_calibration


%%%% Load parameters from calibration
load param_calibrated; 
load eqm_calibrated

%%%% Set up grids for ALPHA and THETA
variable_x          = '\alpha'; range = [0.3,1]; 
ALPHA_grid          = [linspace(0.2,0.4999,20), linspace(0.5, 1, 10)]; 
THETA_grid          = [0.45, 0.50, 0.55];
THETA_grid_legend   = {'0.45','0.50','0.55'};
colors              = [0, 0, 1; 1, 0, 0; 0.9, 0.6, 0.2];
line_style          = {'--';'-';'-.'};

%%%% Options value to be used in fsolve
options                     = optimset('Display','off'); 
options.OptimalityTolerance = 1e-6; 
options.StepTolerance       = 1e-6;
options.MaxIter             = 1000000; 
options.MaxFunEvals         = 1000000; 
options.FunctionTolerance   = 1e-6;

%% 2. --------------------      Now, running comparative statics  ---------
disp('*********************> Comparative Statics <*************************')

par         = par_calibrated;
fun         = define_function_fun(par);     


% model_ should take the following values; 
% 1: Hansen model.
% 2: Hansen model + idiosyncratic productivity + no reallocation
% 3: Hansen model + idiosyncratic productivity + perfect reallocation;
% 4: perfect credit    + no money;
% 5: imperfect credit  + no money;
% 6: imperfect credit  + money (also calibrated parameters)
par_ini = par_calibrated;
do_model_1_and_2; %% Hansen model

for model_ = 3:6 %%%% Not showing Model 3 (close to Model 1)
    %%%% each model file below (e.g., do_model_4) will produce comparative
    %%%% statics results, stored in, e.g., 
    %%%% Consumption_model_4, Inv_model_5, and etc.
    
    switch model_        
        case 3
            disp('-------------> Model 3, idiosyncratic shocks + perfect reallocation<-----------')
            disp('It is the same as Model 1 but with different productivity; therefore we skip.')
            disp(' ')
            
        case 4
            %disp('-------------> Model 4, idiosyncratic shocks + imperfect reallocation + perfect credit <---------')
            par_ini         = par_calibrated;
            par_ini.CHI_q   = 1000; % large number in the credit constraint to represent perfect credit
            do_model_4;
            
            x_grid_model_4 = ALPHA_grid;
            
        case 5
            %disp('-------------> Model 5, idiosyncratic shocks + imperfect reallocation + imperfect credit <---------')         
            par_ini         = par_calibrated;
            par_ini.CHI_k   = 0.20; % this is to distinguish model 5's results better from others     
            do_model_5
            
            x_grid_model_5 = ALPHA_grid;
            
        case 6
            %disp('-------------> Model 6, idiosyncratic shocks + imperfect reallocation + imperfect credit + money <---------')
            par_ini      = par_calibrated;
            fun         = define_function_fun(par_ini);          
            par_ini.CHI_k   = 0.20;
            par_ini.iota    = 0.02;
            plotting_   = 0; % Do not plot Model 6 result individually 
            do_model_6       
    
            x_grid_model_6 = ALPHA_grid;         
    end
end

%% 3. ---------------------  Plot figures
figure()  
subplot(7,3,1)
h = plot(x_grid_model_4, Output_model_4 * 100, 'LineWidth', 2);    
set(h, {'color'}, num2cell(colors, 2), {'LineStyle'}, line_style);
title('Output')
xlabel(variable_x)
ylabel('%')
xlim(range)
ylim([100,115])
grid on

subplot(7,3,4)
h = plot(x_grid_model_4, Inv_model_4 , 'LineWidth', 2);
set(h, {'color'}, num2cell(colors, 2), {'LineStyle'}, line_style);
title('Investment')
xlabel(variable_x)
xlim(range)
ylim([0.17, 0.20])
grid on

subplot(7,3,7)
h= plot(x_grid_model_4, Consumption_model_4, 'LineWidth', 2);
set(h, {'color'}, num2cell(colors, 2), {'LineStyle'}, line_style);
title('Consumption')
xlabel(variable_x)
xlim(range)
ylim([0.57, 0.75])
grid on

subplot(7,3,10)
h = plot(x_grid_model_4, Hours_model_4, 'LineWidth', 2);
set(h, {'color'}, num2cell(colors, 2), {'LineStyle'}, line_style);
title('Employment')
xlabel(variable_x)
xlim(range)
ylim([0.32, 0.34])
grid on

subplot(7,3,13)
h = plot(x_grid_model_4, R_ratio_model_4  * 100, 'LineWidth', 2);
set(h, {'color'}, num2cell(colors, 2), {'LineStyle'}, line_style);
title('R share')
xlabel(variable_x)
ylabel('%')
xlim(range)
ylim([0,100])
grid on

subplot(7,3,16)
h = plot(x_grid_model_4, P_share_model_4  * 100, 'LineWidth', 2);
set(h, {'color'}, num2cell(colors, 2), {'LineStyle'}, line_style);
title('P share')
xlabel(variable_x)
ylabel('%')
xlim(range)
ylim([0,50])
grid on

subplot(7,3,19)
h = plot(x_grid_model_4, Welfare_model_4  * 100, 'LineWidth', 2);
set(h, {'color'}, num2cell(colors, 2), {'LineStyle'}, line_style);
title('Welfare gain (% equiv C)')
xlabel(variable_x)
ylabel('%')
xlim(range)
ylim([-2, 6])
grid on

subplot(7,3,2)
h = plot(x_grid_model_5, Output_model_5 * 100, 'LineWidth', 2);    
set(h, {'color'}, num2cell(colors, 2), {'LineStyle'}, line_style);
title('Output')
xlabel(variable_x)
ylabel('%')
xlim(range)
ylim([100,115])
grid on

subplot(7,3,5)
h = plot(x_grid_model_5, Inv_model_5, 'LineWidth', 2);
set(h, {'color'}, num2cell(colors, 2), {'LineStyle'}, line_style);
title('Investment')
xlabel(variable_x)
xlim(range)
ylim([0.17, 0.20])
grid on

subplot(7,3,8)
h= plot(x_grid_model_5, Consumption_model_5, 'LineWidth', 2);
set(h, {'color'}, num2cell(colors, 2), {'LineStyle'}, line_style);
title('Consumption')
xlabel(variable_x)
xlim(range)
ylim([0.57, 0.75])
grid on

subplot(7,3,11)
h = plot(x_grid_model_5, Hours_model_5, 'LineWidth', 2);
set(h, {'color'}, num2cell(colors, 2), {'LineStyle'}, line_style);
title('Employment')
xlabel(variable_x)
xlim(range)
ylim([0.32, 0.34])
grid on

subplot(7,3,14)
h = plot(x_grid_model_5, R_ratio_model_5  * 100, 'LineWidth', 2);
set(h, {'color'}, num2cell(colors, 2), {'LineStyle'}, line_style);
title('R share')
xlabel(variable_x)
ylabel('%')
xlim(range)
ylim([0,100])
grid on

subplot(7,3,17)
h = plot(x_grid_model_5, P_share_model_5  * 100, 'LineWidth', 2);
set(h, {'color'}, num2cell(colors, 2), {'LineStyle'}, line_style);
title('P share')
xlabel(variable_x)
ylabel('%')
xlim(range)
ylim([0,50])
grid on

subplot(7,3,20)
h = plot(x_grid_model_5, Welfare_model_5  * 100, 'LineWidth', 2);
set(h, {'color'}, num2cell(colors, 2), {'LineStyle'}, line_style);
title('Welfare gain (% equiv C)')
xlabel(variable_x)
ylabel('%')
xlim(range)
ylim([-2, 6])
grid on

subplot(7,3,3)
h = plot(x_grid_model_6, Output_model_6 / min(min(Output_model_6)) * 100, 'LineWidth', 2);    
set(h, {'color'}, num2cell(colors, 2), {'LineStyle'}, line_style);
title('Output')
xlabel(variable_x)
ylabel('%')
xlim(range)
ylim([100, 115])
grid on

subplot(7,3,6)
h = plot(x_grid_model_6, Inv_model_6, 'LineWidth', 2);
set(h, {'color'}, num2cell(colors, 2), {'LineStyle'}, line_style);
title('Investment')
xlabel(variable_x)
xlim(range)
ylim([0.17, 0.20])
grid on

subplot(7,3,9)
h= plot(x_grid_model_6, Consumption_model_6, 'LineWidth', 2);
set(h, {'color'}, num2cell(colors, 2), {'LineStyle'}, line_style);
title('Consumption')
xlabel(variable_x)
xlim(range)
ylim([0.57, 0.75])
grid on

subplot(7,3,12)
h = plot(x_grid_model_6, Hours_model_6, 'LineWidth', 2);
set(h, {'color'}, num2cell(colors, 2), {'LineStyle'}, line_style);
title('Employment')
xlabel(variable_x)
xlim(range)
ylim([0.32, 0.34])
grid on

subplot(7,3,15)
h = plot(x_grid_model_6, R_ratio_model_6  * 100, 'LineWidth', 2);
set(h, {'color'}, num2cell(colors, 2), {'LineStyle'}, line_style);
title('R share')
xlabel(variable_x)
ylabel('%')
xlim(range)
ylim([0,100])
grid on

subplot(7,3,18)
h = plot(x_grid_model_6, P_share_model_6  * 100, 'LineWidth', 2);
set(h, {'color'}, num2cell(colors, 2), {'LineStyle'}, line_style);
title('P share')
xlabel(variable_x)
ylabel('%')
xlim(range)
ylim([0,50])
grid on

subplot(7,3,21)
h = plot(x_grid_model_6, Welfare_model_6  * 100, 'LineWidth', 2); 
set(h, {'color'}, num2cell(colors, 2), {'LineStyle'}, line_style);
title('Welfare gain (% equiv C)')
xlabel(variable_x)
ylabel('%')
xlim(range)
ylim([-2, 6])
grid on

for i = 1:length(THETA_grid)
    entry_legend{i} =  ['\theta = ' THETA_grid_legend{i}];
end
legend(h,entry_legend, 'fontsize',14, 'location','Northeast','Orientation','horizontal');